Libraries:

# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
library(ggordiplots)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------

Source, tranform, select, and relativize MA & LI veg data (hidden)

Species found in at least 5% (11 plots) of plots in Coastal Barrens are: VAPA, GABA, QUIL, GAPR, VAAN, PTAQ, QUCO, PIRI, KAAN, MELI, RUSP, CAPE, QUAL, SMGL, & QUPR

Source and transform envr & categorical data; merge for predictor df to compare nmds against (hidden)

Run NMDS

# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
                   distance = 'bray',
                   k = 3,
                   trymax = 20,
                   autotransform = FALSE)
## Run 0 stress 0.1272962 
## Run 1 stress 0.1275739 
## ... Procrustes: rmse 0.03533478  max resid 0.1292549 
## Run 2 stress 0.1272962 
## ... New best solution
## ... Procrustes: rmse 0.0009770007  max resid 0.003222686 
## ... Similar to previous best
## Run 3 stress 0.1275739 
## ... Procrustes: rmse 0.03479693  max resid 0.1269719 
## Run 4 stress 0.1275743 
## ... Procrustes: rmse 0.03490449  max resid 0.1272927 
## Run 5 stress 0.1272965 
## ... Procrustes: rmse 0.001136387  max resid 0.0037354 
## ... Similar to previous best
## Run 6 stress 0.1275741 
## ... Procrustes: rmse 0.03482774  max resid 0.1271008 
## Run 7 stress 0.1272958 
## ... New best solution
## ... Procrustes: rmse 0.0003339852  max resid 0.00103124 
## ... Similar to previous best
## Run 8 stress 0.1272962 
## ... Procrustes: rmse 0.0006479442  max resid 0.002175514 
## ... Similar to previous best
## Run 9 stress 0.1272961 
## ... Procrustes: rmse 0.0005806935  max resid 0.001964567 
## ... Similar to previous best
## Run 10 stress 0.1272962 
## ... Procrustes: rmse 0.0006353023  max resid 0.002144382 
## ... Similar to previous best
## Run 11 stress 0.1272958 
## ... Procrustes: rmse 0.000133377  max resid 0.0003272971 
## ... Similar to previous best
## Run 12 stress 0.127296 
## ... Procrustes: rmse 0.0005204546  max resid 0.001740775 
## ... Similar to previous best
## Run 13 stress 0.127574 
## ... Procrustes: rmse 0.03500575  max resid 0.1278142 
## Run 14 stress 0.1275736 
## ... Procrustes: rmse 0.03486528  max resid 0.1273871 
## Run 15 stress 0.1272961 
## ... Procrustes: rmse 0.000635704  max resid 0.002127893 
## ... Similar to previous best
## Run 16 stress 0.1275742 
## ... Procrustes: rmse 0.03505738  max resid 0.1279639 
## Run 17 stress 0.1275737 
## ... Procrustes: rmse 0.03491215  max resid 0.1275222 
## Run 18 stress 0.1275735 
## ... Procrustes: rmse 0.0348122  max resid 0.1272177 
## Run 19 stress 0.1275739 
## ... Procrustes: rmse 0.03497939  max resid 0.1277474 
## Run 20 stress 0.127296 
## ... Procrustes: rmse 0.0002188625  max resid 0.0005793351 
## ... Similar to previous best
## *** Best solution repeated 8 times
# converges, stress of 0.127 - 0.128

# 2 dimensions
# mali_2d <- metaMDS(rel_mali5, 
#                    distance = 'bray', 
#                    k = 2,
#                    try = 40,
#                    trymax = 40,
#                    autotransform = FALSE)
# stress still too high to represent in 2 dimensions, lowest is ~ 0.2

Significant species & envr factors for Axis 1 & 2

##            NMDS1       NMDS2 Species  pval
## GAPR -0.43283432 -0.36592469    GAPR 0.026
## VAAN -0.58548007  0.23659024    VAAN 0.004
## QUCO  0.47191797  0.31542190    QUCO 0.010
## PIRI -0.09464334  0.63995157    PIRI 0.001
## KAAN -0.64346682  0.04951010    KAAN 0.001
## MELI -0.21626012 -0.54646216    MELI 0.003
## RUSP -0.62458049  0.06725700    RUSP 0.004
## SMGL  0.06322859  0.74451448    SMGL 0.001
## QUPR -0.58237657 -0.08846107    QUPR 0.011
##                      NMDS1      NMDS2   env.variables  pval
## Time_from_Treat  0.3664852 -0.3459599 Time_from_Treat 0.046
## Moss            -0.1485437  0.5470263            Moss 0.013
## BA_HA            0.5192040 -0.3902518           BA_HA 0.006
## PIRI.BA_HA       0.5939817 -0.4423917      PIRI.BA_HA 0.002

Significant species are: GAPR, VAAN, QUCO, PIRI, KAAN, MELI, RUSP, SMGL, QUPR

Significant environmental factors are: Moss, Basal area per hectare, and PIRI basal area per hectare.

I’m only going to keep one BA measure; PIRI has a lower p value, so thinking that one, but open to other ideas.

3D Plot

This is a 3D ordination being represented by 2D

Using the first two axes, which represent the most variation in the ordination space

GG plot of Axis 1 vs. Axis 2

Plot with significant species & envr factors Axes 1 & 2

Plot with indcator species & envr factors: Axes 1 & 2

Plotting envr factors & species for axes 1 & 3

##             NMDS1      NMDS3 Species  pval
## QUIL  0.001238461  0.6996833    QUIL 0.002
## GAPR -0.413939366  0.3911541    GAPR 0.019
## VAAN -0.571941725 -0.2857426    VAAN 0.003
## QUCO  0.436621926 -0.4459518    QUCO 0.005
## KAAN -0.615704312 -0.2950916    KAAN 0.001
## RUSP -0.572890031 -0.4210261    RUSP 0.001
## QUAL  0.163594982 -0.5134256    QUAL 0.020
## QUPR -0.561433797 -0.2528658    QUPR 0.004
##                NMDS1     NMDS3 env.variables  pval
## avgLD      0.2099211 0.6656505         avgLD 0.003
## BA_HA      0.5201249 0.3122623         BA_HA 0.014
## PIRI.BA_HA 0.6201426 0.2041855    PIRI.BA_HA 0.003

Significant species are: QUIL, GAPR, VAAN, QUCO, KAAN, RUSP, QUAL, QUPR

Significant environmental factors are: average litter depth, Basal area per hectare, and PIRI basal area per hectare.

Graph of Axis 1 vs. Axis 3: indicator species

Graph of Axis 1 vs. Axis 3: significant species

Plotting envr factors & species for Axes 2 & 3

##            NMDS2      NMDS3 Species  pval
## QUIL  0.03601922  0.6998788    QUIL 0.001
## PIRI  0.62635710 -0.2909393    PIRI 0.001
## MELI -0.52040825 -0.2251256    MELI 0.008
## QUAL -0.32982027 -0.5095111    QUAL 0.003
## SMGL  0.74047454  0.1515597    SMGL 0.001
##                   NMDS2      NMDS3 env.variables  pval
## Moss          0.5336551 -0.1538819          Moss 0.018
## Mineral_Soil  0.4287884 -0.3477778  Mineral_Soil 0.013
## Wood          0.1803424  0.4688395          Wood 0.043
## Litter_Duff  -0.2796757 -0.4292763   Litter_Duff 0.041
## avgLD        -0.1887942  0.6465422         avgLD 0.002

Significant species are: QUIL, PIRI, MELI, QUAL, SMGL

Significant environmental factors are: moss, mineral soil, wood, litter-duff, & average litter depth

Graph of Axis 2 vs Axis 3: indicator species

Graph of Axis 2 vs Axis 3: significant species

Run PerMANOVA to look at the impacts of region and treatment on community composition

summary(as.factor(ma.meta.data_veg$Treat_Type))
##  Control   FallRx  Harvest    MowRx SpringRx 
##        5        3        6        3        6
ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)

print(ma.tt)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
##            Df SumOfSqs      R2      F Pr(>F)   
## Treat_Type  4   1.9569 0.31903 2.1082  0.003 **
## Residual   18   4.1769 0.68097                 
## Total      22   6.1339 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise follow up

library(pairwiseAdonis)

pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
## 
## $FallRx_vs_SpringRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.21549 0.10956 0.8613  0.596
## Residual    7  1.75138 0.89044              
## Total       8  1.96687 1.00000              
## 
## $FallRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.48859 0.46035 3.4122    0.1
## Residual    4  0.57276 0.53965              
## Total       5  1.06134 1.00000              
## 
## $FallRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.40379 0.20243 1.7767  0.121
## Residual    7  1.59094 0.79757              
## Total       8  1.99473 1.00000              
## 
## $FallRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.39033 0.23712 1.8649  0.107
## Residual    6  1.25581 0.76288              
## Total       7  1.64614 1.00000              
## 
## $SpringRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.65697 0.28336 2.7678  0.011 *
## Residual    7  1.66151 0.71664                
## Total       8  2.31848 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1   0.4952 0.15598 1.8481  0.042 *
## Residual   10   2.6797 0.84402                
## Total      11   3.1749 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Control
##            Df SumOfSqs     R2      F Pr(>F)
## Treat_Type  1  0.16729 0.0666 0.6422  0.772
## Residual    9  2.34457 0.9334              
## Total      10  2.51186 1.0000              
## 
## $MowRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.72529 0.32578 3.3823  0.017 *
## Residual    7  1.50107 0.67422                
## Total       8  2.22636 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $MowRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.84147 0.41918 4.3303  0.019 *
## Residual    6  1.16594 0.58082                
## Total       7  2.00741 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Harvest_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.58418 0.21102 2.4072  0.038 *
## Residual    9  2.18412 0.78898                
## Total      10  2.76830 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control

Indicator species analysis

library(labdsv)

ma.meta.data_veg$Treat_Group <- NA

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)

mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)

summary(mali_isa)
##      cluster indicator_value probability
## QUPR       2          0.8549       0.003
## KAAN       2          0.7534       0.010
## RUSP       2          0.7097       0.008
## SMGL       4          0.6896       0.012
## 
## Sum of probabilities                 =  4.667 
## 
## Sum of Indicator Values              =  6.68 
## 
## Sum of Significant Indicator Values  =  3.01 
## 
## Number of Significant Indicators     =  4 
## 
## Significant Indicator Distribution
## 
## 2 4 
## 3 1
gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]

print(fridg)
##      group    indval pvalue freq
## QUPR     2 0.8549332  0.003    6
## KAAN     2 0.7534400  0.010   10
## RUSP     2 0.7096551  0.008    7
## SMGL     4 0.6896051  0.012    6
# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is Control

MA_LI MowRx indicators: QUPR, KAAN, RUSP

MA_LI Harvest indicator: SMGL

Finding variation presented by each axis

# axis scores------------------------------
dist1Se <- vegdist(rel_mali5, method = 'bray') # dist in real space
dist2Se <- dist(mali_3d$points[,1], method = "euclidean") # dist. in ordin. space
dist3Se <- dist(mali_3d$points[,2], method = "euclidean")
ax1Se <- lm(dist1Se ~ dist2Se) # distance in real space explained by NMDS axis 1
summary(ax1Se)
## 
## Call:
## lm(formula = dist1Se ~ dist2Se)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31228 -0.06555  0.01538  0.07752  0.28965 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61312    0.01304   47.01   <2e-16 ***
## dist2Se      0.17468    0.01567   11.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1177 on 251 degrees of freedom
## Multiple R-squared:  0.3312, Adjusted R-squared:  0.3286 
## F-statistic: 124.3 on 1 and 251 DF,  p-value: < 2.2e-16
ax2Se <- lm(dist1Se ~ dist2Se + dist3Se) # distance in real space explained by NMDS axis 1 and 2
summary(ax2Se)
## 
## Call:
## lm(formula = dist1Se ~ dist2Se + dist3Se)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.238642 -0.053360  0.004834  0.056850  0.241994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50292    0.01393   36.11   <2e-16 ***
## dist2Se      0.19466    0.01264   15.40   <2e-16 ***
## dist3Se      0.17845    0.01495   11.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09413 on 250 degrees of freedom
## Multiple R-squared:  0.5741, Adjusted R-squared:  0.5707 
## F-statistic: 168.5 on 2 and 250 DF,  p-value: < 2.2e-16
(ax1labSe <- round(summary(ax1Se)$adj.r.squared, 3) * 100) # rsquared of lm of real space explained by axis 1
## [1] 32.9
(ax2labSe <- round(summary(ax2Se)$adj.r.squared * 100 - ax1labSe, 1)) # r squared of lm of real space explaied by axis 1 and 2 less that of axis 1, to get axis 2 alone
## [1] 24.2